Stability of solitary waves in random nonlocal nonlinear media 
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We consider the interplay between nonlocal nonlinearity and randomness for two different nonlin- 
ear Schrodinger models. We show that stability of bright solitons in presence of random perturba- 
tions increases dramatically with the nonlocality-induced finite correlation length of the noise in the 
transverse plane, by means of both numerical simulations and analytical estimates. In fact, solitons 
are practically insensitive to noise when the correlation length of the noise becomes comparable to 
the extent of the wave packet. We characterize soliton stability using two different criteria based 
on the evolution of the Hamiltonian of the soliton and its power. The first criterion allows us to 
estimate a time (or distance) over which the soliton preserves its form. The second criterion gives 
the life-time of the solitary wave packet in terms of its radiative power losses. We derive a simpli- 
fied mean field approach which allows us to calculate the power loss analytically in the physically 
relevant case of weakly correlated noise, which in turn serves as a lower estimate of the life-time for 
correlated noise in general case. 



I. INTRODUCTION 

In diverse physical settings the dynamics of nonlinear 
wave packets is described by the nonlinear Schrodinger 
equation (NLS) including, for instance, nonlinear op- 
tics Q, Bose-Einstein condensates (BEC) [3-5], and wa- 
ter waves 0| . More recently, the NLS equation has been 
also considered in the context of so-called rogue waves [7j . 
The interplay of diffraction/dispersion, which naturally 
tends to spread the wave, and nonlinearity as governed 
by the NLS can lead to formation of solitons - robust 
localized particle-like wave packets, that do not change 
their shape upon temporal evolution or spatial propaga- 
tion n. 

Solitons are ubiquitous in nature and can be found in 
many nonlinear s yste ms ranging from optics Q, physics 
of cold matter [H,[l0| and plasma [U, [l2j to biology [l3| . 
In realistic settings nonlinear systems supporting solitons 
are often subject to random perturbations [14[. Such per- 
turbations may arise from the fluctuation of the external 
linear potential confining the wave, as in case of, e.g., 
BEC's in spat ially and temporally fluctuating trapping 
potentials [15l, KL6j , optical beams in nonlinear dielectric 
waveguides [17j, or waveguide arrays [18j with random 
variation of refractive index, size or waveguide spacing. 
Furthermore, the optical nonlinearity of nematic liquid 
crystals can exhibit stochastic variation due to fluctu- 
ations of the crystal temperature or conditions of sur- 
face anchoring (e.g., roughness) affecting orientation of 
the crystal's molecules pjl, [20|. Similarly, fluctuations of 
temperature will introduce randomness in colloidal sus- 
pensions [2l|, [22| , while noise in the magnetic field em- 
ployed to control scattering length of the Bose-Einstein 
condensate via Feshbach resonance [23| will result in the 
stochasticity of its nonlinear interaction potential [24j . 
Also, stochasticity in the Gross-Pitaevskii equation must 
be often taken into account in order to describe quan- 
tum effects in dilute ultra-cold Bose-gases (e.g. [25]). Fi- 



nally, fabrication-induced spatial fluctuations in periodic 
ferroelectric domain patterns in quadratic media act as 
a source of spatial disorder in the nonlinearity of the 
quasi-phase matched parametric wave interaction affect- 
ing propagation of quadratic solitons [26l428j . 

It has been well appreciated that randomness in the 
linear or nonlinear potential supporting solitons may 
have dramatic consequences on their stability and dy- 
namics depending on the strength of disorder. The pres- 
ence of randomness leads to radiation being emitted by 
the self-guided wave packet, the amount of which de- 
pends crucially on the typical length scale of the fluctu- 
ation or correlation of the noise. The emission of radia- 
tion weakens the self-induced localization and, ultimately 
leads to the decay of solitons. In fact, it has been shown 
that disorder is equivalent to the presence of an effec- 
tive loss in the nonlinear system |27H29j. On the other 
hand, it appears that the interplay between nonlinear- 
ity and weak randomness can lead to diverse interesting 
phenomena, such as random- walk of solitons in the trans- 
verse plane (30l-[32| or Anderson localization [l8|, 0, l33l — 
|35|. Up to now, mostly local nonlinear interaction has 
been considered in studies of solitons in nonlinear ran- 
dom systems. This amounts to Kerr- type nonlinear opti- 
cal response and contact boson interaction in BEC. Re- 
cently however, a few works appeared dealing with ran- 
dom systems that exhibit spatially nonlocal nonlinear- 
ity [HSIiiEI. 

Nonlocality of the nonlinear response appears to be 
common to a great variety of nonlinear systems. Phys- 
ically speaking, nonlocality means that the nonlinear 
response of the medium in a specific location is deter- 
mined by the wave amplitude in a certain neighborhood 
of this location. The extent of this neighborhood is of- 
ten referred to as the degree of nonlocality. Nonlocal- 
ity is common to media where certain transport pro- 
cesses such as heat or charge transfer [37j, diffusion [38| 
and/or drift [39] of atoms are responsible for the non- 
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linearity. It also occurs in media with long-range inter- 
particle interaction. This is the case of nematic liquid 
crystals where nonlinearity involves reorientation of in- 
duced dipoles [So| and in the context of Bose-Einstein 
Condensates with noncontact long-range interatomic in- 
teraction (4ll-[43j . Nonlocality of nonlinearity and its im- 
pact on solitons has been studied extensively in the last 
decade. One of the most important features of nonlocal- 
ity is its ability to arrest catastrophic collapse of mul- 
tidimensional waves I43 - |47| and stabilize complex soli- 
tonic structures |48l-[51|. These stabilizing properties of 
nonlocality have also been identified in the presence of 
randomness. For instance, in recent studies of many- 
soliton interaction in disordered nonlocal medium Conti 
et al have demonstrated that nonlocality leads to for- 
mation of soliton clusters and noise quenching [l9|, [22| . 
Batz et al. [36| reported nonlocality-mediated decrease 
of the quantum phase diffusion and increased coherence 
of quantum solitons while Folli et al. [32| have shown 
that soliton random walk can be efficiently suppressed in 
highly nonlocal media. 

In the present work we will study the effect of non- 
locality on the stability of solitons in nonlinear random 
media. While many previous papers consider only "longi- 
tudinal" random perturbations, i.e., a situation where the 
randomness is only a function of the longitudinal (prop- 
agation) coordinate (e.g. [29]), we will deal here with 
the general case when randomness is a function of both 
propagation and transverse coordinates. We will consider 
two practically relevant models for random nonlocal sys- 
tems. In the first model the randomness contributes ad- 
ditively towards the nonlinear response of the medium. 
In the second model, randomness affects directly parame- 
ters characterizing the nonlinear response of the medium, 
i.e., the noise itself becomes nonlinear. We will show that 
nonlocality stabilizes solitons by effectively increasing the 
correlation length of the random perturbation. By using 
a simplified mean field approach we are able to give a 
lower bound for the life-time of solitary wave packets in 
the case of weakly correlated noise. 

The paper is organized as follows: In Section [III we 
will introduce the afore mentioned nonlocal model with 
additive noise. This model naturally incorporates the 
interplay between randomness and nonlocality, and the 
noise term in the field equation is linear. The different 
effects of randomness and nonlocality on the soliton dy- 
namics will be studied first separately in Subsections III Al 
and IIIB[ and for the full model in III CI In Section ITTT1 we 
investigate our second random nonlocal model, where the 
noise acts multiplicatively and randomness becomes non- 
linear as well. Finally, we will conclude in Section HVl 



II. NONLOCAL MODEL WITH ADDITIVE 
NOISE 

We will consider the evolution of the wave function 
ip(x,t) with x and t denoting generalized transverse 



and longitudinal (propagation) coordinates, respectively. 
The function ijj may represent, e.g., the main electric field 
component of a linearly polarized light beam or the wave 
function of a quantum object such as BEC. We assume 
that ^satisfies the following system of coupled equations 



id t i>(x, t) + d xx ijj(x, t) + pi/j(x, t) = 
-cr 2 d xx p + p = \i 



(la) 
(lb) 



Here, p represents the nonlinear response of the 
medium. In the context of nonlinear optics, p is usually 
identified with a nonlinear refractive index change while 
it accounts for the effective two-body interaction poten- 
tial in case of a BEC. In Eq. (jlb|) , e is assumed to be 
a (5-correlated Langevin noise in both longitudinal and 
transverse coordinates. Then, the noise fulfills (e) = 
and (e(x, t)c(x' , t')) = n 2 5(t - t')8(x - where 



1 N 
lim — fj 



denotes ensemble averaging over different stochastic real- 
izations fj, and n is the so-called coupling strength. By 
solving Eq. (|lb|) via Fourier-transform the system can be 
written as a single nonlocal NLS equation 

id t ^(x,t) + d xx ijj(x,t) 

oo 

+i/>(x,t) / R(x-x') [\ip\ 2 (x',i) + e(x',t)]dx' = 0, 



with the nonlocal response function 

R{x) = ± e -^. 



(2) 



(3) 



Eq. (|lb|) has been used extensively to model the nonlin- 
earity of nematic liquid crystals and soft media including 
colloids [l9|,[32|. Because the noise term e(x, t) is additive 
in the constituent equation, it acts as a source term af- 
fecting the medium independently of whether the actual 
signal (e.g., the optical beam) is present or not. There- 
fore, in the nonlocal NLS Eq. ([2]) noise plays the role 
of a random background potential and is not affected 
by the nonlinearity itself. For instance, in case of ne- 
matic liquid crystals such a noise term may arise due to 
spatial and temporal temperature fluctuations affecting 
linear properties of the medium such as the refractive in- 
dex. In fact, as far as the liquid crystals are concerned, 
random perturbations can also be introduced by fluctua- 
tions of anchoring of the molecules of the crystals at the 
boundaries [20| resulting in random perturbation of the 
molecular orientation across the volume of the crystal. 
The parameter a represents the extent of the nonlocality 
of the nonlinear response and hence defines different non- 
local regimes. Without the noise term, Eq. (pQ) supports 
stable nonlocal solitons 15011 . 
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In Eq. (pQ), the nonlocality parameter a leads to both 
nonlocal nonlinearity and finite correlation length. In 
order to clarify the effect of each of these two constituents 
on stability of solitons we will first discuss both of them 
separately, and then consider their combined action. 



A. Effect of transverse correlation of the noise 

In this subsection we will discuss the effect of trans- 
verse correlation of the randomness on dynamics of local 
solitons. To this end, we will consider the following NLS 



idti/j + d xx ij) + l^lV + v( x , 0^ = 0, 



(4) 



where the random perturbation term is given by the non- 
local relation 



v{x ,t) = jR {x - x >w,tW. 



(5) 



Here, e(x) is again a white noise. For the sake of con- 
sistency with the original model Eq. (pQ), we will use the 
exponential function Eq. (j3]) as kernel function. How- 
ever, we verified that our findings also hold for other, 
e.g. Gaussian, correlation functions. 

The dynamics of localized solutions of Eq. (|4]) is deter- 
mined by two length scales, the width a of the response 
function R and the width 1 / y/\ of the initial soliton so- 
lution 



2A 



cosh y/\i 



(6) 



Figure [T] illustrates typical propagation scenarii of soli- 
tons in this system for successively decreasi ng de gree of 
nonlocality a and coupling-strength n = y/cr/5. Plots 
in the top row [Fig. QJa-b)] represent single stochastic 
realizations for very large a. One can see that in this 
case the peak intensity remains almost constant upon 
propagation. A further decrease of a [Fig. DJc-f)] leads 
to evident decay of the soliton due to emission of radi- 
ation. At the same time, the soliton itself performs a 
random walk in the transverse direction (Gordon-Hauss 
effect [30] ) . The suppression of the soliton random walk 
by nonlocality has been already discussed by Conti et 
al. [32|. Here, we will focus on the noise-induced decay 
of the solitons. 

Since the original noise term e is (5-correlated, one can 
immediately find that 

(r)(x, t)r]{x f ,t f )) = C(x - x')5(t - t'). (7) 

For R given by Eq. (j3j), the function C(x) reads 



C(x) = ^e-M'°(<T+\x\), 



(8) 



and describes the effect of nonlocality on the noise. 
Namely, while the actual random source is represented 
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FIG. 1. Single realizations of the propagation of local solitons 
perturbed by correlated noise for different values of a (degree 
of nonlocality). In (a-f), we have a = 20, 5, 2, 0.5, 0.2, 
0.05, with the respective coupling strength n — y/a/5. The 
decay of the soliton due to the emission of radiation becomes 
stronger for smaller a, even though the coupling strength is 
reduced as well. Apart from the decay, solitons perform a 
random walk in the transverse plane. 



by a white noise the nonlocal character of the nonlinear- 
ity transforms it into an effective colored noise. The non- 
locality acts as a low-pass filter eliminating the high fre- 
quencies of the original noise source and thereby smooth- 
ing out the randomness. This is best appreciated in spa- 
tial Fourier domain (/ = FT[/]) where the initially white 
noise perturbation is modified by a band-pass filter de- 
fined by the Fourier spectrum of the nonlocal response 
function 



fj(k,t) = V27ri?(fc)e(M)- 



(9) 



As a result, the effective noise n exhibits a finite corre- 
lation length determined by the degree of nonlocality a. 
The crucial dependence of the soliton dynamics on a is 
shown in Fig.[T]for single realizations. In the following we 
will investigate the behavior of averaged quantities (•). 



1. random phase-shift 

Let us first consider the situation when the width of 
the nonlocal response function a significantly exceeds the 
spatial extent of the soliton i.e., a ^> 1/a/A, In this case, 



the resulting correlation length of the noise is large, so 
that the dependence of the randomness in the transverse 
plane is basically averaged out, and the noise r] only 
depends on the longitudinal coordinate, rj(x,t) = rj(i). 
Then, the noise term can be effectively removed from 
the original stochastic Eq. (j4]) by the transformation 
ijj(x,t) = i(j(x,t) ex.p[— i f_ rj(t')dt']. In this regime of 
longitudinal-only disorder the soliton maintains its in- 
tensity profile \ip\ 2 while acquiring a random phase-shift. 
In other words, soliton amplitudes will evolve indepen- 
dently of the strength of the noise. However, the random 
phase-shift of the soliton may become apparent in en- 
semble averaged quantities, i.e., decays in time, since 
each realization acquires a different random phase-shift 
upon propagation. 

Generally, when the noise rj(x,t) depends on both 
transverse and longitudinal coordinates x and t, it is pos- 
sible to show that (see App. [A] for details) 



id t + d xx + i 



C(0) 



W + <MV> = 0. (10) 



Unfortunately, due to the nonlinear term (|^| 2 ?/;) it is 
not possible to solve Eq. ([TO]) directly. However, for large 
correlation length a the perturbation to the soliton due to 
spatial noise is weak and one can approximate d xx (^) + 
(\ifj \ 2 ifj) w A where A is the soliton parameter from 
Eq. ©. Then, Eq. (TH) gives 



2 e -C(0)t 



(11) 



for times t sufficiently small. Physically speaking, in 
Eq. (pTJ) we neglect random walk and radiative decay 
of the soliton. 

To verify the above findings we must resort to numer- 
ical analysis. In our simulations of Eq. (J2J) we use a 
semi-implicit method in the interaction picture [52| and 
absorbing boundary conditions to solve the stochastic 
partial differential equation. The random term is imple- 
mented by using Gaussian distributed random numbers 
generated by a Box-Muller method [53| . Due to the noise- 
induced emission of radiation the absorbing boundaries 
emulate parts of the wave function ip that leave the fi- 
nite numerical box and thus lead to an effective decrease 
of power (or number of particles) P = J |^| 2 dx in our 
simulations. P would remain a conserved quantity if we 
could use an infinitely large numerical box. 

Throughout this subsection we will consider ipo(x) = 
V2/cosh(x) (A = 1) as an initial condition. Fig. [2fa) 
depicts the evolution of \(ip(x = 0,t))| 2 (solid lines), for 
different values of the correlation length a. The auto- 
correlation C(0) = n 2 /4a is fixed to 1/20 by adjusting 
the coupling strength n, so that there is only one predic- 
tion (dashed line) for the evolution \(ip(x = 0))| 2 from 
Eq. ([11]). As expected, we find that Eq. JTTJ is exact for 
infinite correlation length a = oo (i.e., no spatial noise). 
For finite values of <r, Eq. ([TT]) holds for small times t 
only, because of the combined effect of emission of ra- 
diation and random walk. In fact, for large a > 1 it is 
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FIG. 2. (a) Evolution of \{i/j(x = 0)),£| 2 (solid lines) com- 
pared to the analytical prediction Eq. (|lip (dashed line). The 
correlation length is a = oo, 3, 1, and 1/128 for the green, 
blue, black and red curve, respectively; the auto-correlation 
C(0) = n 2 /4cr = 1/20 is kept constant, (b) The correspond- 
ing soliton random walk (solid lines) compared to the analyt- 
ical prediction of Eq. (|12[) (dashed lines). In our simulations, 
ensemble averages are taken over 128 realizations. 



mainly the random walk which causes the deviation of 
\{ijj{x = 0))| 2 from the predictions of Eq. ([TT]) , whereas 
for small a < 1 radiative losses to the soliton are more 
important (see also Fig. [T]). 



2. random walk 



The random walk of the solitons can be estimated an- 
alytically. With X(t) = fx\ip\ 2 dx/P, SX(t) = X- (X), 
it is possible to show that [32| 



(5X(tf) 



16nV 
~3Pg 



-t 6 



(12a) 



« = fjj M*i) [dxM*i)] R ( x i - x s) (12b) 
x ^ (x 2 ) [^^0(^2)] R(x2 ~ x 3 )dxidx 2 dx 3 . 



In Fig. [2fb) we compare predictions of Eq. ([T2]) (dashed 
lines) to numerical results (solid lines). Whereas we find 
perfect agreement for large degree of nonlocality a, there 
is a clear deviation for smaller a. This is an effect of 
the increasing radiation losses with smaller a. When 
the soliton power decreases upon propagation, the ran- 
dom walk decreases as well. To illustrate this behav- 
ior, let us consider the situation when a is significantly 
smaller than the width of the spatial soliton un- 

der consideration. We see that according to Eq. (fT2j) 
the strength of the random- walk only depends on n 2 , be- 
cause k/Pq — » A 3//2 /15 for a — y 0. However, we have to 
take into account that in this regime radiation losses to 
the soliton affect k/Pq during propagation and introduce 
deviations from (SX(t) 2 ) oc t 3 p54f. 
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3. Hamiltonian 

Finding an analytical estimate for the deformation of 
the soliton due to randomness is difficult. However, the 
stabilizing effect of correlation can be illustrated by con- 
sidering the quantity 



(13) 



i.e., the Hamiltonian of Eq. (4) when r] = 0. Thus, 
in the limit n — >• the Hamiltonian % is a conserved 
quantity, but it becomes time dependent for finite cou- 
pling strength. Several papers [55l457| already empha- 
sized the fact that (H) is then a linear function of time t, 
(H(t)) = Ho + jt. In fact, it is possible to compute the 
ascent 7 analytically (see App. IB]). Then, the final result 
for the time evolution of the averaged Hamiltonian reads 



(H) = Ho 



d 2 C{x) 
dx 2 



Pt = H 



x=0 



n 2 P 
4a 3 1 



(14) 



We assume that the soliton does not change its shape sig- 
nificantly until the critical time t-u = — H0/7 is reached, 
where (H(t<u)) = and 7 = n 2 P/4a 3 . This formula 
already emphasizes the huge impact of the nonlocality 
mediated correlation length a on the soliton propaga- 
tion: the critical time tu scales with a 3 , whereas the 
noise amplitude enters only as 1/n 2 [see also Fig. [3f b)] . 
In Fig. 03^a) we compare predictions of Eq. (fT4|) (dashed 
lines) to numerical results (solid lines). As long as the 
total power P is a conserved quantity in the simulations, 
i.e., the numerical box is large enough, we find perfect 
agreement. However, noise generated radiation spreads 
quite fast and will eventually reach any numerical bound- 
ary. Moreover, the radiation gives important contribu- 
tions to the Hamiltonian %, because the gradient norm 
J \ Vip\ 2 dx is very sensitive to small-scale fluctuations. 
Even though we used numerical boxes of size 300-600 
and 8000-16000 points in x, for larger propagation times 
t numerical curves deviate from the ideal linear behavior. 



4- radiative losses 

As we have seen so far, the effect of emission of ra- 
diation and subsequent loss of power is crucial for the 
soliton dynamics, particularly when a is small compared 
to 1/vX In the following, we will refer to soliton sta- 
bility in terms of decay or decrease of the peak intensity 
l^lmaxW = max cc |^l 2 upon propagation. Due to the 
stochastic character of the propagation Eq. (jl]) we will 
discuss, in fact, the ensemble averaged quantity ( 1^1 max)- 

It is actually possible to estimate analytically the ra- 
diative losses to the soliton in the limit a <C 1/vX To 
this end, let us assume that the noise term rjip in Eq. (j4]) 
acts perturbatively (of the order S <C 1) on the soliton 
^0, and ip = (-00 + x) ex P(^A£), where x(#,£) is of the 
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FIG. 3. (a) Evolution of Hamiltonian H [Eq. (H3l) ] 

(solid lines) compared to analytical predictions from Eq. (|14p 
(dashed lines). The correlation length is a = 3, 1, and 0.5 
for the purple, black and red curve, respectively; the auto- 
correlation C(0) = n 2 /4a — 1/20 is kept constant, (b) shows 
the analytical prediction for tn (dashed black line) as a func- 
tion of a. Numerical values are indicated by circles, where 
the colored circles correspond to the colors used in (a). In 
our simulations, ensemble averages are taken over 128 real- 
izations. 



order S <C 1 as well. Then, in order S 1 we find 

idtX - Ax + d xxX + 2|^o| 2 X + V^oX* + #0 = 0. 

At initial time t = we have the pure soliton ^0 without 
any perturbation, and thus x( x > t = 0) = 0. For small 
times At we can therefore write down a formal solution 
for the perturbation 



At 



x(x,At)&i / rj(x,t)i/;o(x)dt. 



(15) 



When the correlation length a of the noise is much 
smaller than the width of the soliton 1/V% the pertur- 
bation x is spectrally much broader than the soliton ip. 
Therefore, we can assume that x w iU essentially describe 
radiation, or, in other words, the part of the total wave 
function which is completely alien to the soliton. In or- 
der to proceed, we will compute the ensemble average of 
P^ f = J \x(x, At)| 2 dx, i.e., the power of the radiation 
produced in the time interval At: 



At At 



r)(x, t)^o(x)r](x, t')2pQ(x)dt'dtdx} 







= C(0) J |^o(^)| 2 dxAt = C(0)P At, 

where we used Eq. flTJ). On the other hand, Eq. (jl]) is con- 
servative which dictates that the fraction of power con- 
verted to radiation P^ is lost to the soliton. It is known 
that the Schrodinger soliton can adapt adiabatically to 
losses by moving along the family branch towards lower 
powers [58|, i.e., the soliton parameter \(t) and power 
Po(t) oc y/\(t) become slowly decreasing in time. If we 
assume that the radiation, once produced, does not in- 
teract anymore with the soliton and disperses quickly, 
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FIG. 4. (a) Time t 1 / e where the soliton power drops be- 
low Po(t — 0)/e versus degree of nonlocality a /as- Numer- 
ical simulations confirm the analytical prediction Eq. (|17[) , 
ti/e — m the weakly nonlocal limit a /as <C 1. The 

dotted line represents an eye-guide oc (a/as) 3 ^ 2 . The auto- 
correlation (7(0) = n 2 /4a — 1/20 is kept constant, (b) shows 
the corresponding evolution of the power of the solitons. The 
correlation length is<j = oo, 3, 1, 0.5, and 1/128 for the violet, 
black, green and red curve, respectively. In our simulations, 
ensemble averages are taken over 128 realizations. 



Eq. ([T5]) becomes valid for any interval [£, t + At], and we 
can conclude that 



(Po(t + At)) = (P o (t))[l-C(0)At]. 



(16) 



Here, Po(t) ex y/X(i) denotes the power in the solitonic 
part of the total wave function ip(x,t). With At — >• we 
therefore get 

(Po(t)) = P (t = 0)e~ c ^\ (X(t)) = X(t = 0)e- 2C ^ t . 

(17) 

Interestingly, we find that noise induced radiation losses 
to the soliton are described by the same constant (7(0) 
we already found in the context of the "random phase- 
shift" in Eq. ([TU]) . However, the two effects are en- 
tirely different: Eq. ([TP]) describes | (ip(x, t)) | 2 , whereas 
Eq. (pTTj) approximates the soliton power related to 
(\ip(x — X(t),£)| 2 ), i.e., when both random phase-shift 
and random walk play no role. 

Figure 3] illustrates the decay of soliton power Pq due 
to radiation. In Figure H^a), we show the time t\i e for 
the ensemble averaged soliton power integrated over an 
interval of 5as centered around the center of the soliton 
X(t). Here, as denotes the initial 1/e width of the soliton 
intensity profile, i.e., as = 1.085/vX Figure |4jb) shows 
the dynamics of the power decay for different values of 
the correlation length a, and confirms that our simplified 
model Eq. (fT7j) works well for sufficiently small a. Com- 
paring the times ty obtained in the previous subsection 
and ti/ e , we see that ty <C t\j e . Depending on the con- 
text and interest, one has to decide which one to use. 
Clearly, if ty is large the soliton will propagate for long 
times without significant change. The time t\i e is useful 
when we are interested in estimating the destruction of 
the wave packet, i.e., ask the question how long will the 
soliton survive. 



We note that using a similar reasoning is possible to es- 
timate radiation losses for an arbitrary initial condition 
ip(x,t = 0), provided that it is spectrally narrow com- 
pared to the noise and random walk is negligible. Then, 
it is possible to write down an evolution equation for the 
so-called mean field ^mf = (^(x,t)), 



id* 



, / ,2 .C(0) 



MF 



o, 



(18) 



and we have I^mf 2 ~ (\ip(x,t)\ 2 ). We checked 
the validity of Eq. (jTHJ) numerically for various non- 
solitonic initial data iMx,t = 0). Using perturbative 
variational techniques J59[ , it is simple to show that soli- 
tons in Eq. (fT8|) decay with the same rate as that found 
in Eq. (IT7|) . 



5. white noise (a ^ 0) 

Before moving over towards nonlocal solitons, let us 
discuss the limit of purely white spatial noise, i.e., when 
a — >• 0. In this case the spatial noise term in Eq. (pQ) has 
to be interpreted appropriately, depending on the under- 
lying physics. In principle, as a represents a relevant 
physical length scale (60| , one has to make sure that it is 
resolved numerically by choosing a sufficiently fine mesh 
with spatial step-size Ax <C a. However, this may re- 
quire high computational costs for small a. For example, 
in numerical simulations employing a = 1/128 we have 
to use at least Ax « 10 -3 . In order to reduce compu- 
tational costs, one could say that this noise is effectively 
^-correlated. Then, according to Eq. flSJ, in the limit 
a — >> and fixed coupling strength n (i.e., the strength of 
the random walk is kept constant) radiation losses C(0) 
in Eq. ([TT]) formally go to infinity. From a physical point 
of view this is not a problem because a may be small 
but never zero. However, in the case that we do not (or 
can not) resolve a, e.g., in a macroscopic approach, the 
noise becomes formally ^-correlated, R(x) = 5(x) and 
C{x) =n 2 S(x). 

In our numerical scheme (5-correlation means that 
(7(0) = n 2 / Ax. Provided that the wave function ip is suf- 
ficiently resolved on the mesh with step-size Ax ^> a, nu- 
merical simulations using ^-correlated noise are supposed 
to mimic those employing the actual response R(x). Fig- 
ure [5] shows simulation results for ^-correlated noise. Pre- 
dictions obtained from Eq. (fT7|) [or, alternatively, from 
the mean field Eq. ([15]) ] are in excellent agreement with 
the numerical results. Thus, as far as radiation losses are 
concerned, in order to mimic an actual response R(x) 
with finite a, we have to choose the same value for the 
auto-correlation (7(0). This essentially means that we 
have to use an effective coupling strength n e ff ex V Ax 
for the (5-correlated noise. Indeed, comparing the curves 
for a = 1/128 in Fig. 2] with those in Fig. [5] where 
n = V0.05Ax (i.e., (7(0) = 1/20), we see that it is actu- 
ally possible to substitute noise with extremely small cor- 
relation length by ^-correlated noise and thus obtain the 



7 




FIG. 5. (a) Evolution of the ensemble averaged peak inten- 
sity ( 1^1 max) (solid red lines) compared to analytical predic- 
tions Eq. (fT7|) (dashed black lines) for ^-correlated white noise. 
T he coup ling strength is n = 0.005, 0.01, 0.02, 0.03, and 
V0.05Ax, respectively; the spatial step size Ax = 80/2048 is 
kept constant. The last value n — V0.05Ax is chosen such 
that it yields the same decay rate (7(0) for the soliton power as 
the case a = 1/128 in Fig.[H Graph (b) shows the correspond- 
ing evolution of the power of the solitons. In our simulations, 
ensemble averages are taken over 128 realizations. 



same results with much less computational costs. How- 
ever, by doing so we sacrifice the accurate description of 
the random walk of the solitons: For a coupling strength 
n e ff oc \J Ax random walk becomes mesh dependent and 
vanishes for Ax — >> [see Eq. (fl2]) ]. 

Obviously, here we have to make a choice which effect 
we want to describe correctly and independently of the 
step-size Ax in the J-correlated case. On one hand, if 
we choose the random walk to be grid independent (n 2 
fixed), (7(0) and therefore radiation losses become grid 
dependent and go to infinity for Ax — >• 0. This strategy 
was used, e.g., in Ref. [61]. On the other hand, if we want 
to describe radiation losses correctly, we have to fix the 
auto-correlation (7(0) to the value of its counterpart for 
the original correlated noise. 



B. Effect of nonlocality of nonlinear potential 

The stabilizing effect of nonlocal nonlinearities with 
respect to collapse (44l-[47j. perturbations of the initial 
soliton profil e I45I , and even support of higher-order soli- 
tons [48145 ll [62| has been extensively discussed in the 
literature. Based on these results one might be tempted 
to expect that the nonlocal character of the self-induced 
nonlinear potential itself will be sufficient to weaken the 
destabilizing effect of (5-correlated random perturbations 
on the wave function by spatially smoothing out the ran- 
dom perturbation. To clarify this issue we will analyze 
in this subsection the following nonlocal model 



idtil* + d xx ^ + [R * H 2 ] ^ + ei/> = 0, 



(19) 



i.e., opposite to the model considered in the previous sub- 
section, the nonlocal response function R gets convoluted 
(indicated by *) with the nonlinearity \t/j\ 2 and not with 



the white noise e(x,t). Thus, the nonlocality of the sys- 
tem will affect only the nonlinear self-induced potential 
without modifying the noise term which will remain 5- 
correlated. 

As in the previous subsection, it turns out that again 
the evolution of the average peak intensity of the soliton 
(l^l 2 ) can be accurately described by a mean field quan- 



tity |^; 



MF 



J ), and the mean field amplitude ^mf 



is governed by the following equation 



id t + d xx + R* |^mf| 



; C(0) 



^MF 0, 



(20) 



in complete analogy to Eq. (fTHj) . This finding indicates 
that loss of power (or particle number) and subsequent 
decay of the effect of self-trapping of the soliton solely 
depends on (7(0) = v? j Ax. In particular, for given noise 
level all solitons exhibit the same rate of power loss in- 
dependently of the degree of nonlocality of the nonlin- 
ear self-trapping potential, which is somehow counter- 
intuitive since in many other situations a stabilizing effect 
of nonlocal nonlinear potentials has been observed. 

Detailed numerical investigations (not shown) confirm 
that nonlocality of the nonlinear response as represented 
in Eq. (fT9]l has no effect on the decay of solitons in the 
presence of white noise. In fact, the only difference be- 
tween local and nonlocal nonlinear response is, in this 
respect, the modification of the soliton transverse pro- 
file. 



C. Full nonlocal model 

In the complete model Eq. (pQ) nonlocality affects both 
the self-induced nonlinear potential and the spectrum of 
the noise. In the last subsections, these two aspects have 
been considered separately. We will now discuss the soli- 
ton dynamics following from the full nonlocal model of 
nonlinear media Eq. (pQ). In the following consideration 
we fix the initial power of the soliton Po(t = 0) = 4. 
Then, the width of the initial soliton as is just a func- 
tion of the nonlocal length a, which is depicted in the 
inset of Fig. [6jc). The time-evolution of the analogous 
Hamiltonian for different values of a 



1 



(21) 



is shown in Fig. (6]^a). It turns out that the analytical re- 
sults obtained in Subsection III Al hold, i.e., with Furutsu- 
Donsker-Novikov we find again Eq. (|T4|) . The only noted 
difference is that the initial Hamiltonian Ho (as well as 
the soliton shape function) is now dependent on a, and 
has to be calculated by using the numerical exact non- 
local soliton profile. As a consequence, the dashed lines 
representing the analytical predictions in Fig. [6fa) have 
different starting points at t = (solid lines show their 
numerical verifications). Moreover, the resulting "life- 
time" tq-i = —41-Loo~ 3 /n 2 P is slightly smaller than in the 
previous case [see Fig. EJc)]. 
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FIG. 6. Nonlinear nonlocal model with additive noise Eq. JT]). 

(a) Evolution of the Hamiltonian % [Eq. (|2T|) ] (solid lines) 
compared to analytical predictions from Eq. ([T4|) (dashed 
lines). The nonlocal length is a = 3, 1, 0.5, and 1/128 for the 
purple, black, green, and red curve, respectively; the auto- 
correlation C(0) — n 2 /4a = 1/20 is kept constant, (b) shows 
the corresponding evolution of the power of the solitons. (c) 
shows the analytical prediction for tn (dashed black line) as a 
function of a. Numerical values are indicated by circles, where 
the colored circles correspond to the colors used in (a) and 

(b) . The inset in (c) shows the functional dependency of the 
relative correlation length a /as on the nonlocal length a for 
Po = 4. (d) shows the time t 1 / e where the soliton power drops 
below Po(t = 0)/e versus degree of nonlocality a /as- Numer- 
ical simulations confirm the analytical prediction Eq. (|22[) , 
ti/e — 1/C(0), in the weakly nonlocal limit a /as <C 1. The 
dotted line represents an eye-guide oc (a/as) 2 . In our simu- 
lations, ensemble averages are taken over 128 realizations. 



Finally, let us consider the effect of radiation losses 
to the nonlocal solitons in Eq. (pQ). Figure EJb) shows 
the ensemble averaged soliton power, integrated over an 
interval of 5as centered around the barry center X(t), 
for different values of a. In the weakly nonlocal limit, 
we recover the results found in the previous subsections 
[compare to Fig. |4^b) and Fig. Efb)], i.e., an exponential 
decay of the ensemble averaged soliton power 



<P (t)> = P (f = 0)e 



-C(0)t 



(22) 



with decay rate C(0) [cf. Eq. (fTTj)]. FigureEJd) illustrates 
the resulting 1/e life-time of the soliton as a function of 
the degree of nonlocality <j/<js- We also verified that for 
spectrally narrow arbitrary initial conditions ip(x,t = 0) 
the mean field equation ([20]) holds. 



III. NONLOCAL MODEL WITH 
MULTIPLICATIVE NOISE 

In this section we will briefly discuss the interplay of 
randomness and nonlocality in a second nonlocal model. 
We consider a nonlinear medium described by the follow- 
ing coupled system 



id t il>(x, i) + d xx ip(x, t) + pip(x, i) = 



(23a) 



-(J 2 d xx p + p=(l + 6)\^\ 2 . (23b) 

Unlike the previously discussed model with additive noise 
here the noise amplitude is a function of the soliton in- 
tensity. As before, e(x, t) is assumed to be a J-correlated 
Langevin noise in both longitudinal and transverse coor- 
dinates. Such a model may represent an inhomogeneous 
nonlinearity of liquid crystals, colloids or a nonlinear- 
ity arising from imperfections in periodic poling in quasi 
phase matched nonlinear structures. The above coupled 
system is equivalent to the following nonlocal Schrodinger 
equation 

idt^J + d xx i) + R(x) * [(1 + e(x, t))\^j\ 2 }^ = 0, (24) 

with the usual exponential kernel R. 

In the multiplicative model Eq. (|23]h the effective 
strength of the noise depends crucially on the power of 
the soliton, since the noise term is multiplied by the in- 
tensity. Similarly to the case of additive noise [Eq. (PQ)], 
we can estimate radiative losses in the weakly nonlocal 
limit: Analogously to Eq. ([T6]h one finds that the power 
of the soliton is governed by 



(P (t + At)) = (P (*)> 



n 2 At 

4a 



t)\ 6 dx. (25) 



Using the soliton profile, Eq. (j6j) [63|, we find the evolu- 
tion equation for the soliton parameter 



cVA 
~dT 



8r^ 



(26) 



which can be solved analytically with the initial condition 
A(0) = 1, 



\(t) = 



15 



A/480n 2 t/o- + 225' 



(27) 



The corresponding expression for the soliton power can 
be found by P = 4\/A. 

Figure EJb) shows the evolution of ensemble averaged 
soliton power obtained from numerical simulations. In 
the case of spectrally broad noise, i.e., small a, we find 
perfect agreement with our analytical results Eq. (f27j) . 
In order to be able to compare with previous findings, 
n 2 /4<7 = 1/20 has been fixed to the same value. The 
main difference between additive and multiplicative noise 
is that in the latter case the radiative decay of the soli- 
tons is not exponential. Therefore, we consider the time 
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^75%? where the ensemble averaged soliton power drops 
below 75% of Po(t = 0). As expected, Fig.[7Jd) reveals a 
significant increase of t 75 % with increasing nonlocality a. 

Unlike in the case of additive noise, we cannot write 
down a mean field equation or give an analytic expression 
for the Hamiltonian H [Eq. (|2T|) ]. since the effective noise 
depends now on ip itself in a nontrivial way. The times 
t-u •> where (H) becomes zero in our simulations are de- 
picted in Fig. (He). Also the random- walk now depends 
on the evolution of the soliton. Assuming naively, that 
the soliton profile is conserved upon propagation, so that 
we can take the initial profile as spatial noise filter for all 
times, one can find that the random walk of solitons is 
described by a formula analogous to Eq. (fT2j) 



(SX(t) 2 ) 



16n 2 K 
~3P! 



-t a 



(28a) 



^o(zi) [d x ^o(xi)\ R(xi - x 3 ) (28b) 
x i/jo(x 2 ) [d x i/jo(x 2 )} R(x 2 - x 3 )|^o(^3)| 4 cl^icl^2dx 3 



In case of a large nonlocal length, Eq. ([28]) yields good 
results, whereas it fails for smaller a due to significant 
changes in the profile and radiation upon propagation 
[see Fig. [7(a)]. For better readability of the figure, we do 
not show the weak random walk of the case a = 1/128 
in Fig. [7(a). 



IV. CONCLUSIONS 

In this paper we investigated systematically the inter- 
play between nonlocality and randomness for two differ- 
ent prototypical physical models. Our first model consid- 
ers the case, where random fluctuations are present inde- 
pendently of the nonlinear wave (additive noise), whereas 
in the second model randomness itself is dependent on 
the wave function (multiplicative noise). Even though 
we restrict our analysis to an exponential nonlocal func- 
tion, our results should hold for arbitrary kernels. We 
were particularly interested in stability of solitons. The 
stochasticity in the constituent equation leads to radia- 
tion which causes the soliton to loose power (or norm). 
In order to characterize the time-scales of this radiative 
power loss, we derived two criteria. The first one involves 
the Hamiltonian of the noiseless part of the model equa- 
tions, and gives and estimate for how long the soliton can 
propagate without changing its shape. It turns out that 
this "life-time" of the soliton depends crucially on the 
degree of nonlocality, i.e., the ratio of correlation length 
of the noise over width of the soliton. Even though we fix 
the auto-correlation of the noise when changing the de- 
gree of nonlocality, which means that the noise amplitude 
increases for larger nonlocal length a, we find a dramatic 
enhancement of the soliton life-time through nonlocal- 
ity. Interestingly, however, the main impact on soliton 
stability comes from the nonlocality induced increase of 
the correlation length and not from the nonlocality of the 
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FIG. 7. Nonlinear nonlocal model with multiplicative noise 
Eq. (|23[) . (a) Soliton random walk (solid lines) compared to 
the analytical prediction of Eq. ([28]) (dashed lines), (b) shows 
the corresponding evolution of the power of the solitons. The 
correlation length is a = oo, 3, 1, and 1/128 for the green, 
blue, black and red curve, respectively; the coupling strength 
n is adjusted such that n 2 /4a = 1/20 is kept constant, (c) 
shows numerical values of tn as a function of a. Colored 
circles correspond to the colors used in (a) and (b). The 
inset in (c) shows the functional dependency of the relative 
correlation length a /as on the nonlocal length a for Po = 4. 
(d) shows numerical values of the time £75% where the soliton 
power drops below 0.75Po(£ = 0) versus degree of nonlocality 
a /as- In our simulations, ensemble averages are taken over 
128 realizations. 



nonlinearity. Our second criterion for soliton stability in- 
corporates the fact that a soliton prone to radiative losses 
can adiabatically transform itself into another member of 
its family with lower power. As a result, the wave packet 
stays confined upon propagation, but gradually decreases 
its guided power. Therefore, it makes sense to discuss, 
e.g., a "1/e-life-time" in terms of power. It turns out 
that the "1/e-life-time" as well strongly increases with 
the nonlocal length <j. In case of spectrally very broad 
noise, resp. very small correlation length compared to 
the extent of the wave packet, we found analytical ex- 
pressions quantitatively describing the loss of power due 
to radiation in both nonlocal models. These analytical 
results are particularly useful derive upper estimates for 
decay rates of the solitons. In particular, the convergence 
of weakly correlated noise to formally ^-correlated noise 
has been discussed in depth. Finally, we addressed spa- 
tial random- walk of the solitons modified by radiation as 
an additional important dynamical feature. 
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Appendix A: Derivation of the averaged equation 

For the derivation of the ensemble averaged Eq. (|T0l) we 
resort to the Furutsu-Donsker-Novikov formula |l4l l64l- 



I (^^){v(x,tHx',t'))dx'dt' (Al) 



Causality in time is reflected by the upper integration 
limit. This formula involves the variational derivative of 
the wave function if; with respect to the noise term 77. To 
compute this quantity we write Eq. (jl]) in integral form 



+ i J [d xx + \i>(x,t')\ 2 +r 1 (x,t')]iP(x,t')dt', 



(A2) 



where ipo(x) is the initial condition at t = 0. Then, with 
^0|y = S{x - x')5(t - t') (see, e.g., Q) we find for 
< t' < t 



5ip(x, t) 



5r)(x',t') 



y^j = i J [d xx + rj(x,t 



Sr}(x',t') 



+ i 



I 



Sr](x , ,t / ) 



dt" + ii/;(x,t f )5(x-x f ). 



Here we used the causality principle again, namely that 
5 ^f2 = for t" < t f . Finally, in the limit t f -> t we 
get 



(A3) 



(A4) 



and Eq. ([TP]) can be found in straight forward manner by 

taking the ensemble average of Eq. dU . 

Appendix B: Derivation of the averaged Hamiltonian 

The ascent 7 of the ensemble averaged Hamiltonian 
(H) equals the time-derivative of the (H). Thus, we have 
to compute 



d t (H) = d t (l \d^\ 2 - -M 4 dx) 



J 



(irj {ijj*d xx il) - ijjd xx ijj*))&x 



Here, * means complex conjugation, and for the time- 
derivatives of the wave function ^(x,i) we plugged in 
Eq. (jlj. In the next step, we make use of the Furutsu- 
Donsker-Novikov formula (IA1I): 



t 00 00 



d t (H)=i 



,5[il>*(x,t)d xx fl>(x,t)] 



-OO —OO —OO 



Srj(x' 



Sr ] (x , ,t / ) ' 



-UI< 



S[i>*d xx il> - ipd xx <ip* 
Sr](x / ,t) 



■}C(x-x')dxdx' 



In the second step, we performed the time-integration 
over t 1 . By evaluating the variational derivatives Eq. (|A3[) 
and integration by parts we read 

dt(n) = ±C(0)( J tl>*d xx *l> + c.c. dx) 

--([[ *P*5(x - x')d xx [C(x - x')^} + c.c. dx'dx). 



Most of the integrals above turn out to be zero, which 
can again be seen by integration by parts, and we find 



d t (H) 



d 2 C(x) 
dx 2 



([ \4>(x,t)\ 2 dx). 
=0 J 
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